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ABSTRACT 

We present the result of iV-body simulations of dynamical evolution of triple 
massive blackhole (BH) systems in galactic nuclei. We found that in most cases 
two of the three BHs merge through gravitational wave (GW) radiation in the 
timescale much shorter than the Hubble time, before ejecting one BH through a 
slingshot. In order for a binary BH to merge before ejecting out the third one, it 
has to become highly eccentric since the gravitational wave timescale would be 
much longer than the Hubble time unless the eccentricity is very high. We found 
that two mechanisms drive the increase of the eccentricity of the binary. One 
is the strong binary-single BH interaction resulting in the thermalization of the 
eccentricity. The second is the Kozai mechanism which drives the cyclic change of 
the inclination and eccentricity of the inner binary of a stable hierarchical triple 
system. Our result implies that many of supermassive blackholes are binaries. 

Subject headings: black hole: physics - - black hole: binary — galaxies: nu- 
clei — stellar dynamics — Three-body problem:general — Gravitational Wave 
Radiation : LISA — methods: n-body simulations 

1. Introduction 

When two galaxies merge, the central BHs sink toward the center of the merger remnant 
and form a binary system (Begelman, Blandford, & Rees 1980, henceforth BBR). Whether or 
not this binary BH can merge in the timescale shorter than the Hubble time is an important 
question for the understanding of the evolution and growth of BHs and galactic nuclei. 
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BBR predicted that the hardening of the BH binary through the stellar dynamical 
interaction with field stars (FSs) would slow down after the loss cone is depleted. After the 
depletion, the growth will be driven by the refilling of the loss cone by two-body relaxation, 
and a simple theoretical estimate results in the merger timescale orders of magnitude larger 
than the Hubble time. 

Whether or not this slowing down of evolution actually occur has been studied with 
iV-body simulations (Quinlan & Hernquist 1997; Milosavljevic & Merritt 2003), Until very 
recently, however, there has been serious discrepancy between theory and simulation results, 
and also among simulation results. Theory predicts that the timescale must be proportional 
to the number of particles used to express the parent galaxy, since the relaxation timescale 
is (neglecting the weak logiV term) proportional to the number of particles in the system 
N. The results of iV-body simulations range from no dependence on iV (Chatterjee et al. 
2003; Milosavljevic & Merritt 2001) to iV 1/3 (Makino 1997). 

Makino & Funato (2004) performed iV-body simulations with much larger number of 
particles (up to 10 6 ) and for longer time compared to previous works, and concluded that 
the iV-dependence of the evolution timescale of BH binary is consistent with the theoretical 
prediction. Berczik, Merritt, & Spurzem (2005) obtained similar result, by performing even 
longer simulations for somewhat smaller iV (up to 256k). 

Thus, as far as the stellar dynamical effects are concerned, BH binary, once depleted 
the loss cone, would not merge within the Hubble time, unless there remain lots of gas in 
the central region of merger remnants. 

If BHs in galaxies, in particular giant ellipticals, are mostly binaries, growth of the BH 
through merging might be difficult. If one galaxy with binary BH and the other with single 
BH merge, the central BHs form a triple system, which is generally unstable. The most 
likely outcome would be the ejection of at least one BH, preventing the growth of the central 
BH binary. In fact, triple interaction of BHs was first considered as an mechanism to eject 
BHs from the center of a galaxy(Saslaw et al. 1974). The timescale for the merging through 
GW radiation is given by 



where mg is the mass of BHs (assuming equal masses), a and e are the semi-major axis 
and eccentricity of the binary. In the case of a BH with 10 8 M Q , the orbital separation must 




(1) 



9(e) 



(1 - e 2 ) 7/2 



(2) 



1 + (73/24)e 2 + (37/96)e' 



,4' 



- 3- 



become less than 0.1 pc in order for the merging timescale to become shorter than the Hubble 
time, if we assume the average thermalized value of e = 0.7. 

Since the gravitational wave timescale depends on the eccentricity very strongly, if some 
mechanism drives the binary BH to high eccentricity the merging timescale might become 
short. Makino & Ebisuzaki (1994) argued that the maximum eccentricity a binary BH would 
experience is much higher than the average value, since a binary BH must undergo around 
10 strong single-binary scattering events before ejecting out a single BH. They concluded, 
using simple analytic arguments, that the binary BH in a triple BH system would merge, 
for the case of typical BH of 10 8 M Q in a typical elliptical galaxy within the Hubble time. 
Another possibility is the change of eccentricity through Kozai mechanism(Kozai 1962). 

In this paper, we investigate the evolution of triple BH systems by means of direct TV- 
body simulations. Our initial model consists of an iV-body realization of a spherical galaxy 
and three BH particles. We started with hierarchical configuration in which two BH particles 
initially form a binary. The third BH particle is placed at some distance to the binary. We 
changed the initial orbit of the third particle, to see its effect on the evolution of the binary. 
We incorporated the effect of GW radiation on the relative orbits of BH particles. 

We found that if the initial orbits of three BHs are coplanar and if the total angular 
momentum is small, the binary BH hardens through repeated binary-single BH interactions. 
The interaction stops either when one BH is ejected out of the galaxy or when the BH binary 
merges through GW radiation. Which of the two dominates depends on the depth of the 
potential well of the parent galaxy. Merging is more likely for deeper potential. For typical 
giant ellipticals, potential is deep enough that merging dominates. 

When the orbits are not coplanar, hierarchical triples are formed and the inner binary 
evolves through Kozai mechanism, in a fair fraction of cases. Kozai mechanism leads the 
large-amplitude oscillation of eccentricity. Thus, merging is more likely in three-dimensional 
configuration than in coplanar systems. 

The structure of this paper is as follows. In section 2, we describe the initial models 
and method of our numerical simulation. In section 3, we show the result of the cases in 
which three black holes are initially in coplanar orbits. In section 4, we show the result of 
more general cases of non-coplanar initial conditions. We summarize our result of numerical 
simulation in section 5. 
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2. Initial Models and Numerical Methods 



2.1. Models 



The initial setup is schematically shown in figure 1. We placed three BH particles in 
iV-body models of a spherical galaxy. The quantitative properties of models and initial 
conditions of our simulation is summarized in table 1. 

For the galaxy model, we used a King model with Wo = 7, where Wo is the nondi- 
mensional central potential of King models (King 1966; Binney & Tremaine 1987). We 
performed several runs changing the number of particles in a galaxy from 16K to 128K. 
We adopted the the standard iV-body units (Heggie & Mathieu 1986), in which M m i = 
1.0, E ga i = — 0.25, G = 1. Here, M m i and E ga i are the total mass and total binding energy 
of the galaxy and G is the gravitational constant. The relation to the physical units will be 
discussed in section 2.5. 

We placed three blackholes in the parent galaxy model. All blackhole particles have the 
same mass of 0.01 in the iV-body units. Two blackholes are initially placed as a binary with 
eccentricity e ~ 0.7 in the center of the galaxy model. The orbital plane of this binary is the 
x — y plane, the third blackhole is placed at position (1,0,0). 

We have changed the initial velocity of the third particles in several different ways. In 
the simplest case, it is placed at rest. In one series, we gave initial velocity along y axis, and 
in another series along z and in the third one in y — z plane. Tables 2 and 3 list all runs 
we performed. We scaled the positions and velocities of stars so that the total energy of the 
system, including that of black hole particles, becomes —1/4. 



2.2. Equations of Motion 

The equation of motion for field particles is given by 

d 2 r f:i 



dt 2 



F //,i + F fB,H ( 3 ) 



where F// ;i and FfB,i are acceleration due to field stars and BH particles respectively. They 
are given by 
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Fig. 1. — Initial setup of BH particles in the parent galaxy is shown schematically. Small 
dots and large filled circles denote field particles and BH particles, respectively. The arrow 
shows the initial velocity of the 3rd blackhole. i is the inclination. 
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Here m/^ and r^j are the mass and the position of field star with index i, and rriB,i, rs,i 
are those of black hole particle i. We used softening parameters e// and e/B for field star- 
field star interactions and field star-black hole interactions. Their values are 6ff = 0.01 and 
efB — 10~ 6 , respectively. 

The equation of motion for BH particles is 



where Fsf,i and Fbb,i are acceleration due to field stars and BH particles, respectively. They 
are given by 



We set the softening between BH particles as €bb = 0. The second term, F GW , is the 
quadrupole approximation for the orbital change due to GW radiation (Damour 1987), which 
is expressed as 



dt 2 



— Ffi/,t + F B B,i, 
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v — v BB,ij — v B,i — v B,j, 



r — r BB,ij — r B,i — r B,j, 



(12) 
(13) 



where Vb^ is the velocity of black hole i. 
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2.3. Numerical Method 

We integrated the equations of motions using the 4th order Hermite scheme (Makino & 
Aarseth 1992) with individual variable time step. In order to apply the Hermite scheme, we 
need the time derivative of the acceleration. For the gravitational wave term (Eq. 11), we 
used the following form as the time derivative: 




a — a BB,ij — a B,i ~ a B,j, (15) 

where ae,i is ith black hole acceleration. 

In order to calculate the acceleration due to field particles, we used GRAPE-6 (Makino 
et al. 2003), the special purpose computer for the gravitational A-body problem. Forces 
from BH particles, both Newtonian and gravitational wave terms, are calculated on the host 
computer. 

In all runs, the energy (corrected for the loss of energy through GW radiation)is con- 
served to 0.1%. 

2.4. Correspondent Physical Scales 

We interpret the mass unit in our simulation as 10 10 M Q . Thus, the mass of BH particles 
is 10 8 M Q . The total mass of 1O 1O M corresponds to the mass of the nucleus of a galaxy 
containing a BH. In other words, the galaxy model we use is not really the model of the 
entire galaxy, but only that of the central region. The reason we limit our model to the 
central region is simply to reduce the number of particles A, in other words, to reduce the 
calculation cost. This modification means that the potential depth of our model galaxy is 
somewhat less than that of a real galaxy with the same velocity dispersion. To see the effect 
of this difference, we have changed the velocity dispersion of the galaxy and investigated 
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the effect. In the standard runs, we set the velocity dispersion of the initial King model 
as 300 km/s. In other words, we set the light velocity c in the iV-body unit to 706. We 
tried two different values of velocity dispersion, 600 km/s and 150 km/s, or the value of the 
light velocity 353 and 1412, respectively. Note that by changing the value of c actually we 
change the strength of the GW term in equation (11) and thus change the evolution of the 
BH particles. 

3. Result I — Simplest configuration 

In this section, we describe the results of runs in which the third BH is initially at 
rest, to illustrate the typical behavior and the dependence of the results on the artificial 
parameters of numerical experiments such as the softening and the number of particles. We 
discuss the results of more general configurations in section 4. 

Figure 2 shows the evolution of the semi-major axis a and the eccentricity e of the BH 
binary in one run of model 64kW7 (hereafter called run 64kW7A). In this and all other 
figures, the semi-major axis and eccentricity are defined as those of the most strongly bound 
pair of BH particles. Initially (for time T < 1.8), a shows smooth decrease. This is due 
to the dynamical friction from field stars. The jumps in a and e at T ~ 1.8 is the result 
of first strong encounter between the binary and the third body. In this case, the binary 
becomes wider and more eccentric. The three BH particles experienced quite complex series 
of interactions until T ~ 10. At this point, the third BH was ejected out, and it took some 
time before it came back and interacted strongly with the binary at T ~ 13. After this 
interaction, the eccentricity of the binary reached 0.991, resulting in the quick orbital decay 
by GW emission. 

In this case, the high eccentricity is realized only after the last strong interaction. How- 
ever, there are several other periods when the eccentricity was large, like around T = 7 
(e ~ 0.97) and T = 9 (e ~ 0.94). High values of eccentricity are not very unusual after 
strong interaction events. 

In figure 3, orbits of three blackholes in run 64kW7A projected on the x — y plane 
are shown. The top panel shows the trajectory of two BH particles which are initially in a 
binary, for the period of < T < 1.8. The third BH particle is out of the plotting region and 
thus not shown. Initially, the orbits decay quickly. As a result the orbital period becomes 
shorter and the change in one orbit becomes smaller. In the meantime, the center of mass 
of the binary moves along the y axis, because of the gravitational force from the third BH 
particle which was initially located at (0,1,0) and falling down to the center of the galaxy. 
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Fig. 3. — Trajectories of the BHs before three body encounter(T = — 1.8) and just before 
and after the last strong encounter (T = 12.85 — 12.95). 
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The bottom panel of figure 3 shows the last strong encounter. The binary BH (of BHs 
numbered 1 and 3) comes in from the right side, and the third BH (BH2) comes in from the 
left side. After the complex encounter, a new binary of BHs 1 and 2, with the very high 
eccentricity, flies out to upper-right corner and BH3 to the lower-left corner. 

3.1. Merging Probability 

The last column of table 2 shows the fraction of runs in which merger occurred. The 
divisor is the total number of runs performed for the model and dividend is the number of runs 
in which merging occured. For runs of model 64kW7 with velocity dispersion a = 300km/s, 
in three out of five runs the BH binary merged before single BH is kicked out of the parent 
galaxy. In two other cases, a single BH was kicked out. For runs of model 64kW7High 
(a = 600km/s), two out of three resulted in merging. All of runs of model 64kW7Low 
(a = 150km/s) resulted in ejection. 

The total potential depth of typical giant ellipticals is somewhere between galaxy models 
64kW7 and 64kW7High. Thus, we can conclude that a triple BH system, if formed in giant 
ellipticals, is likely to end up in merging of two of the three BHs, leaving out a binary BH 
system. THe Ejection of one BH is not impossible, though. In the remaining part of this 
section, we check the effect of artificial parameters such as the softening and the number of 
particles. 

3.2. Effect of the softening 

In figure 4, we show the evolution of semi-major axis of binaries of two runs with 
different softening parameters between FS-FS. Solid and dashed curves correspond to runs 
with ejf = 0.01 and 0.001, respectively. Though the details are different, the overall evolution 
is similar. We performed five runs with e// = 0.001, and one runs ended up in merging. This 
result is consistent with that for e// = 0.01. 

3.3. Effect of the number of particles 

In figure 5, we show the evolution of semi-major axis of binaries in runs with different 
number of particles for the parent galaxy. Again, the overall behavior is not much different. 
For N = 16k and N = 128k, we again performed ten and five runs, and four runs each ended 
up in merging. We can conclude N has little effect on the result. 
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Table 1. Model Parameters 



Parameter 


Symbol 


Value 


Mass of galaxy 


M-gal 


1 n 
l.U 


Mass of FS 


rrif 


1/16384- 1/131072 


Mass of BH 


m B 


0.01 


Number of FS 


N f 


16- 128/c 


Number of BH 


N B 


3 


Gravitational constant 


G 


1 


Total energy 


Egal 


-0.25 


Softening between BH and BH 




0.0 


Softening between BH and FS 


e B f 


10~ 6 


Softening between FS and FS 


e ff 


0.01-0.001 




2 4 6 8 10 12 14 

T 

Fig. 4. — The evolution of semi-major axis for runs with different values of the softening 
parameter eff for field particles. Solid and dashed curves are the results of runs with eff = 
0.01 and 0.001, respectively. All other parameters are the same as those in model 64kW7. 
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Fig. 5. — The evolution of semi-major axis for different values of N. Solid, long-dashed and 
short-dashed curves are the results of runs with N = 16k, 64k and 128k, respectively. All 
other parameters are the same as those in model 64kW7. 
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This result is quite different from that for the evolution of two massive BHs (one BH 
binary) in the central region of the galaxy, where the number of particles, i.e., the relaxation 
timescale, determines the evolution timescale of the BH binary In the case of two BHs, the 
loss cone is soon depleted, and the only way for the field particles to interact with the BH 
binary is to enter the loss cone through two-body relaxation. Thus, it is quite natural that 
the evolution timescale depends on N. However, in the case of a triple BH system, three BHs 
alone determines the evolution, and field particles works just as the background potential 
and the source for the dynamical friction. Thus, there is nothing like loss-cone depletion. 



3.4. Effect of the initial galaxy model 

In addition to models with Wo = 7, we tried three runs with Wo = 11. As can be 
seen in table 2, there is not much difference in the final outcome. This result looks a bit 
counter-intuitive, since, as we have seen in section 3.1, the merging probability does depend 
on the depth of the potential relative to the speed of the light. 

The reason why the outcome is not much different in cases with Wq = 7 and 11 is that 
initial binary BH hardens much more quickly in models with Wq = 11 than in models with 
W = 7. As a result, ejection of BH occurs after smaller number of triple interactions and 
the chance to reach high eccentricity is smaller for higher central density. This difference 
more than compensates the deeper potential. 



4. Result II — General configurations 

In section 3, we considered the simplest case in which a binary is at the center of the 
galaxy and the third BH particle falls from around the half-mass radius. This is clearly a 



Table 2. Models for section 3 



Model name 


N f 


Wo 


e ff 




c 


merging probability 


16kW7 


16k 


7 


0.01 


300km/s 


706 


4/10 


32kW7 


32k 


7 


0.01 


300km/s 


706 


2/5 


64kW7 


64k 


7 


0.01 


300km/s 


706 


3/5 


128kW7 


128k 


7 


0.01 


300km/s 


706 


4/5 


64kW7High 


64k 


7 


0.01 


600km/s 


353 


2/3 


64kW7Low 


64k 


7 


0.01 


150km/s 


1412 


0/3 


64kWll 


64k 


11 


0.001 


300km/s 


677 


1/3 


64kW7Eps 


64k 


7 


0.001 


300km/s 


706 


1/5 
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rather special case, and the result may be biased. In this section we consider several different 
configurations. 

Figure 6 shows the result of runs in which the third BH has the initial velocity perpen- 
dicular to the orbital plane of the binary BH (runs 64kI90Vx, see Table 3). The case with 
v z = 0.1 is qualitatively similar to the result of free-fall runs (figure 2). Orbital elements 
of the binary shows noncontinuous changes after strong encounters with the third BH. The 
time interval for the encounter seems to be longer, due to the larger total angular momentum 
of the triple BH system. 

The case with v z = 0.3 shows quite different behavior. The eccentricity shows periodic- 
like change for T = 11 ~ 18, but the semi-major axis does not show much change during the 
same period. This behavior is driven by the Kozai mechanism. The period and magnitude 
of the oscillation change in time, because the orbital elements of both the inner and outer 
binaries changes due to the interaction with other field stars. Around T = 18, the hierarchical 
system of three BHs became unstable, and strong three-body interactions followed. 

In the case of v z — 0.7, the inner binary merged in the first cycle of the Kozai cycles. 

Figure 7 shows the results of three runs with different initial inclination of the orbital 
plane of the third BH. As we've seen in figure 6, in the extreme case of large v z perpendicular 
to the binary orbital plane, the end result is the merging in the first Kozai cycle. In the case 
of % — 7r/3, Kozai cycle drives oscillations of large amplitude, and semi- major axis shrinks 
each time when the eccentricity reached the peak. In the case of % = tt/6, BHs neither merge 
nor escape until T = 200. Kozai cycles are also observed. The maximum eccentricity in one 
cycle is, however, not high enough to drive the merging of inner binary. In the case of the 
coplanar orbits, Kozai cycle of a small magnitude is observed for T = 20 to 90. Since the 
initial orbit is coplanar, there cannot be any Kozai cycle if three BHs remain exactly in the 
same orbital plane. However, since the number of particles is finite, fluctuation in the forces 
from field particles introduces random changes in the inclination of inner and outer binaries, 
resulting in the configuration for which the Kozai cycle can show up. 

To summarize the results of 3D triple hierarchical triple BH system is often 

formed when the third BH falls in. In such a triple system, a stable oscillation in eccentricity 
(and inclination) driven by Kozai mechanism is often observed. In a fair fraction of the cases, 
this oscillation results in the merging of the inner binary by GW emission. 
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Fig. 6. — The evolution of the semi-major axis a and eccentricity e for runs with non-zero 
initial V z of the third BH particle. Left top, right top, and bottom panels show the results 
for V z = 0.1,0.3,0.7, respectively. 
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o.i 

0.03 
0.01 

0.003 
0.001 



0.003 



i=60 inner 
i=60 outer 




0.1 
0.03 
0.01 
0.003 
0.001 



i=30 inner 
i=30 outer 



0.3 



0.1 



0.03 



0.01 



10 20 30 40 50 60 70 80 



i=30 inner 

i=30 outer ----- 

20 40 60 80 100 120 140 160 180 200 




0.003 



Fig. 7. — The evolution of the semi-major axis a and eccentricity e for runs with initial 
\V\ = 0.7 for the third BH particle. Initial velocity of the third BH is in yz plane and the 
angle with the xy plane is 60 (left top), 30 (right top) and (bottom) degree. Each panel 
shows the semi-major axis and the eccentricity of inner (thick solid line) and outer binaries 
(thin dashed line). 



5. Discussion 



5.1. Escaping cases and free-floating BHs 

We have seen that the merging of the BH binary is more likely than the ejection of 
single BH (or also the binary BH). However, the ejection is not a rare event. Out of the five 
runs with N = 64k and o = 300km/s listed in table 2, two runs ended up in the ejection. 
Though this result is somewhat biased by the shallow galactic potential, runs with deeper 
potential still gives a fair number of escaping events. 

Thus, our result suggests that, though in most cases the binary BH do merge, there 
may be cases in which one or more BHs are kicked out. In other words, there may be a fair 
number of massive BHs free-floating in intergalactic field. 

An interesting question is if such a floating BH can be observed. One possibility is 
through strong lensing, but since the Einstein radius, when placed around z = 0.5, is around 
10~ 2 arcsec, it would be hard to detect it even with HST. Another possibility is the detection 
of X-ray emission by Bondi accretion of intracluster hot gas, if BH is in a cluster of galaxy. 
The luminosity L is given by 

L = riiiiBC 2 , (16) 

where rj is the radiative efficiency, is the accretion rate and c is the light velocity. The 
accretion rate is given by 

rh B = F(7)7rG' 2 m|p 00 /a^ , (17) 

where ^(7) is a function of adiabatic index of the gas 7, and are density and sound 
speed at infinity, respectively (Bondi 1952). We take rj = 0.1 and ^(7) = 1. For a BH of 
1O 8 M in the Virgo cluster, whose typical hot gas density and sound speed is given by 

Tl 

»- ~ 17 x 10 1.0x10- w^ cm - (18) 

kT \V2 

4.3 xlO 7 ' 



2Sv) Cm/S ' (19) 



V2keV 

Table 3. Models for section 4 



Model name 


N f 


W 


i 


\v\ 


a 


c 


final state 


64kI90Vl 


64k 


7 


90 


0.1 


300km/s 


706 


escape 


64kI90V3 


64k 


7 


90 


0.3 


300km/s 


706 


escape 


64kI90V7 


64k 


7 


90 


0.7 


300km/s 


703 


merge 


64kI60V7 


64k 


7 


60 


0.7 


300km/s 


703 


merge 


64M30V7 


64k 


7 


30 


0.7 


300km/s 


703 


neither merge nor escape by T = 200 


64kI0V7 


64k 


7 





0.7 


300km/s 


703 


merge 
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where n and kT are the number density and temperature of the intracluster hot gas and we 
adopted n — 1.0 x 10~ 4 cm~ 3 and kT = 2keV (Kikuchi et al. 2000) as typical values, thus 
the luminosity would be 

L mrgo = 9.0 x 10 19 m B erg/s (20) 

Thus, it can be observed as a moderately luminosity X-ray source. The question is how it 
can be discriminated from AGNs. 



5.2. Conclusion 

In this paper, we investigated the dynamical evolution of triple black hole systems in the 
central region of a galaxy. We found that in most of cases two black holes merge through GW 
radiation. The merging timescale by GW radiation is greatly reduced by high eccentricity of 
the binary BH. The high eccentricity is realized by two different mechanisms. First one is the 
random change of the eccentricity after binary-single BH interaction as suggested by Makino 
& Ebisuzaki (1994), and the other is the Kozai mechanism(Kozai 1962). Even so, there 
may be a fair number of free-floating BHs ejected out of the parent galaxy by a slingshot. 
They might be observed as X-ray sources, if they stay within a cluster of galaxies. 

We thank Tatsushi Matsubayashi and Toshiyuki Fukushige for stimulating discussions 
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library. 
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